home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / ellipse.i < prev    next >
Text File  |  1995-07-26  |  2KB  |  86 lines
  1. /*
  2.    ELLIPSE.I
  3.    Complete elliptic integrals.
  4.    (Maybe someday can find cn, sn, dn algorithms...)
  5.  
  6.    $Id$
  7.  */
  8. /*    Copyright (c) 1994.  The Regents of the University of California.
  9.                     All rights reserved.  */
  10.  
  11. /* ------------------------------------------------------------------------ */
  12.  
  13. func EllipticK(k)
  14. /* DOCUMENT EllipticK(k)
  15.    return E(k), the complete Elliptic Function 
  16.    integral from 0 to pi/2 of  1/sqrt(1-k^2 *(sin(x))^2) dx
  17.    uses the arithmetic/geometric mean method. (Abramowitz & Stegun p598)
  18.    k must lie in   -1 < k <1.
  19.  
  20.    E(k) diverges as log(4/sqrt(1-k^2)) as k->1
  21.  
  22.    SEE ALSO: EllipticE
  23.  */
  24. {
  25.   EPS= 1.E-15;
  26.   a= 1.0;
  27.   b= sqrt(1.0-k*k);
  28.   c= abs(k);
  29.   do {
  30.     an= 0.5*(a+b);
  31.     bn= sqrt(a*b);
  32.     cn= 0.5*(a-b);
  33.     a= an;
  34.     b= bn;
  35.     c= cn;
  36.   } while (anyof(abs(c)>EPS));
  37.  
  38.   return 0.5*pi/a;
  39. }
  40.  
  41. func EllipticE(k)
  42. /* DOCUMENT EllipticE(k)
  43.      return E(k), the complete Elliptic Function
  44.      integral from 0 to pi/2 of  sqrt(1-k^2 *(sin(x))^2) dx
  45.      k must lie in -1<= k <=1
  46.  
  47.      Uses the arithmetic/geometric mean method. (Abramowitz & Stegun p598) 
  48.  
  49.    SEE ALSO: EllipticK
  50.  */
  51. {
  52.   EPS= 1.E-15;
  53.  
  54.   /* Separate out k=1 since K diverges and d->1 */
  55.   mask= (abs(k)==1.0);
  56.  
  57.   list= where(mask);
  58.   if (numberof(list)) E1= array(1.0,numberof(list));
  59.  
  60.   list= where(!mask);
  61.   if (numberof(list)) {
  62.     c= k(list);
  63.     a= 1.0;
  64.     b= sqrt(1.0-c*c);
  65.     d= 0.0;
  66.     twofact= 0.5;
  67.     do {
  68.       an= 0.5*(a+b);
  69.       bn= sqrt(a*b);
  70.       cn= 0.5*(a-b);
  71.       dn= d + twofact*c*c;
  72.       twofact*= 2.0;
  73.       a= an;
  74.       b= bn;
  75.       c= cn;
  76.       d= dn;
  77.     } while (anyof(abs(c)>EPS));
  78.     En1= 0.5*pi*(1.0-d)/a;
  79.   }
  80.  
  81.   /* merge k=1 and k!=1 results */
  82.   return merge(E1,En1,mask);
  83. }
  84.  
  85. /* ------------------------------------------------------------------------ */
  86.